#ifndef __mml_primitives_h__
#define __mml_primitives_h__

#include <iostream>

namespace mml
{
  /** Maximum element. */
  template <typename Type>
  Type const& vector_max(int n, Type const* x);

  /** Minimum element. */
  template <typename Type>
  Type const& vector_min(int n, Type const* x);

  /** Vector cross product. */
  template <typename Type>
  void vector_cross(Type const* a, Type const* b, Type* c);

  /** Print a vector. */
  template <typename Type>
  void vector_print(int n, Type const* data, std::ostream& s = std::cout);

  /** Read a vector. */
  template <typename Type>
  void vector_read(int n, Type* data, std::istream& s);

  /** Copy a vector */
  template <typename Type>
  void vector_copy(int n, Type* target, Type const* source);

  /** The infinity norm (magnitude of largest component). */
  template <typename Type>
  Type vector_norminf(int n, Type const* x);

  /** The squared 2-norm of a vector. */
  template <typename Type>
  Type vector_normsq(int n, Type const* x);

  /** The 2-norm of a vector. */
  template <typename Type>
  Type vector_norm(int n, Type const* x);

  /** Angle between two vectors. */
  template <typename Type>
  Type vector_angle(int n, Type const* a, Type const* b);

  /** Normal to a vector. */
  template <typename Type>
  void vector_normal(int n, Type const* x, Type* y);

  /** Interchange contents of two (non-overlapping) vectors. */
  template <typename Type>
  void vector_swap(int n, Type* x, Type* y);

  /** Normalize a vector. */
  template <typename Type>
  void vector_normalize(int n, Type* x);

  /** Inner product of two vectors. */
  template <typename Type>
  Type vector_ip(int n, Type const* x, Type const* y);

  /** Scale a vector. */
  template <typename Type>
  void vector_scale(int n, Type* x, Type const& scale);

  /** Add a scaled vector to another. */
  template <typename Type>
  void vector_add(int n, Type* x, Type const* y, Type const& scale);

  /** Add a vector to another. */
  template <typename Type>
  void vector_add(int n, Type* x, Type const* y);

  /** Subtract a vector from another. */
  template <typename Type>
  void vector_subtract(int n, Type* x, Type const* y);

  /** Negate a vector. */
  template <typename Type>
  void vector_negate(int n, Type* x);

  /** Set a vector to a value. */
  template <typename Type>
  void vector_set(int n, Type* x, Type const& s);

  /** Reset a vector to default value. */
  template <typename Type>
  void vector_reset(int n, Type* x);

  // Calculate the difference ||a||^2 * ||b||^2 - (a * b)^2 in a numerically stable way.
  template <typename Type>
  Type vector_abmabsq(int n, Type const* a, Type const* b);

  // Calculate the difference ||a|| ||b|| - a * b in a numerically stable way
  template <typename Type>
  Type vector_abmab(int n, Type const* a, Type const* b);

  /** Vector outer product. */
  template <typename Type>
  void vector_outerProduct(int m, int n, Type const* a, Type const* b, Type* A);



  /** Grassmann product of two quaternions. */
  template <typename Type>
  void quat_grassmann(Type const* a, Type const* b, Type* c);



  /** Print a matrix. */
  template <typename Type>
  void matrix_print(int m, int n, Type const* A, std::ostream& s);

  /** Generate the transpose of a matrix. */
  template <typename Type>
  void matrix_transpose(int m, int n, Type const* A, Type* At);

  /** Multiply a vector by a matrix b = Ax. */
  template <typename Type>
  void matrix_matvecmul(int m, int n, Type const* A, Type const* x, Type* b);

  /** Multiply a matrix by a matrix C = A * B. */
  template <typename Type>
  void matrix_matmatmul(int m, int n, int o, Type const* A, Type const* B, Type* C);

  /** Reset a matrix. */
  template <typename Type>
  void matrix_reset(int m, int n, Type* matrix);

  /** Add a matrix to another. */
  template <typename Type>
  void matrix_add(int m, int n, Type* target, Type const* source);
}


#endif
